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We propose a model for co-evolving ecosystems that takes into account two levels of description 
of an organism, for instance genotype and phenotype. Performance at the macroscopic level forces 
mutations at the microscopic level. These, in turn, affect the dynamics of the macroscopic variables. 
In some regions of parameter space, the system self-organises into a state with localised activity and 
power law distributions. 
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Complex extended systems showing a lack of scale in 
their features appear to be widespread in nature, being 
as diverse as earthquakes , creep phenomena |^ , ma- 
terial fracturing fluid displacement in porous me- 
dia 0,^, interface growthj? 
and biological evolution [[14 



Tc| ] , river networks 1 11 - ^] 
ITf . At variance with cqui- 
hbrium statistical mechanics, these systems do not need 
any fine tuning of a parameter to be in a critical state. 
To explain this behaviour, Bak, Tang and Wiesenfeld in- 
troduced the concept of self-organised criticality (SOC) 
through the simple sand-pile model [|l8|,^. In recent 
years, several models with extremal dynamics have been 
shown to exhibit SOC in the presence of a sufflciently 
decorrelated signal pOf . 

Several years ago, Bak and Sneppen (BS) proposed a 
SOC model for the co-evolution of species. In the BS 
model each species occupies a site on a lattice. Each site 
j is assigned a fitness, namely a random number between 
and 1. At each time step in the simulation the smallest 
fitness is found. Then the fitnesses of the minimum and 
of its nearest neighbours are updated according to the 
rule 
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that assigns a new fitness /n+i at time rt -I- 1 to the cho- 
sen lattice site. This corresponds to the extinction of the 
less fit and its impact on the ecosystem. Indeed, in the 
original BS model, the function F is just a random func- 
tion with a uniform distribution between and 1. The 
system reaches a stationary critical state in which the 
distribution of fitnesses is zero below a certain threshold 
and uniform above it (the actual value of the threshold 
depends on the updating rule). It has been shown [ p2|j23| 
that the exact nature of the updating rule is not relevant. 
Indeed, the use of a chaotic map instead of a random up- 
date, preserves the universality class (even if the final 
distribution may be altered). 

As a result of the dynamical rules, the BS system 
exhibits sequences of causally connected evolutionary 
events called avalanches [Ell . The number of avalanches 



N follows a power law distribution N{s) s where s is 
the size of the avalanche and r ~ 1.07 [|^|4|j2|]. Other 
quantities also exhibit a power law behaviour with their 
own critical exponents. Prominent among them are the 
first and all return time distributions of activity (a site 
is defined as active when its fitness is the minimum one) , 
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where t/ - 1.58 and Ta - 0.42 |2C|,|2^,g5|| . 

In the BS model, two basic ingredients are needed for 
SOC to occur [||: 

i) Order from extremal dynamics (minimum rule). 

ii) Disorder from the updating rule (stochastic or other- 
wise). 

An oversimplification in the BS model is apparent: 
each species is described by a single variable. The mini- 
mum rule is applied to this variable, and so is the effect of 
mutations. In natural evolving systems, however, at least 
two — interacting — levels of organisation are present, and 
both play a role in the evolution. Mutations occur at a 
microscopic level, namely at the molecular level of the 
genotype. This affects a macroscopic level, defining the 
phenotype. Natural selection acts on the phenotype, al- 
lowing the survival of some and the extinction of others, 
according to the observed power law distributions. 

In this paper we propose a new model, namely an Ex- 
tremal Coupled Map Lattice (ECML), that takes into 
account this two level structure. The ecosystem is rep- 
resented as an ensemble of N species arranged on a one- 
dimensional lattice. Each species is described by means 
of macroscopic variable Xi subject to a nonlinear dynam- 
ics and coupled to its nearest neighbours. In general, 
Xi can be identified with a population, or some other 
function of the phenotype. The control parameter of 
the nonlinear dynamics is identified with the microscopic 
level (genotype). To find the exact function connecting 
these two levels of description is beyond the scope of this 
simplified model. For this reason, we chose a logistic map 
for the independent evolution of Xi with Xi acting as the 
nonlinear parameter in the map [|2q| . With this in mind. 
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the evolution of each species is given by: 

<+i = (1 - + I [fi<- 
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where f{x) is the logistic map f{x) = Xx{l — x) (in gen- 
eral any chaotic function should do). Each site has its 
own Ai, extracted from a fixed distribution g(A) [ p7| . 
The local coupling of strength a emulates the ecologi- 
cal interaction between neighbouring species. A general 
description should include inhomogeneities in e. Here, 
for simplicity we limit our analysis to the homogeneous 
case, ei = e. 

We consider this CML as the substrate on which evolu- 
tion takes place. The parameter Ai, that determines the 
behaviour of Xi, is regarded as the microscopic level, the 
genotype, and is subject to mutation. We suppose that, 
through mutation, a species is able to alter this strategy 
to adapt to the environment defined by the collective be- 
haviour. For this we propose an extremal mechanism, 
akin to the BS model ||21]. The species that, at each 
time step, has the minimum value of a;, is considered the 
candidate to mutation. Its A is replaced by a new value 
drawn from the distribution g{X). 

We can summarise these simple dynamical rules by 

1. Evaluate expression (^). 

2. Find the site with the absolute minimum fitness on 
the lattice (this site will be called the "active" site) . 

3. Change the value of A of the active site. 

4. Go to step 1. 




FIG. 1. Space-time diagram of the active sites in a CML 
of 500 sites. Each dot is an active site, i.e. the site with 
the minimum x at a certain time. AA — 0.1 in all four 
cases. (a)Ao = 3.7, e — 0.5. (b)Ao — 3.9, e = 0.5. (c) Ao — 3.8, 
e = 0.1. (d)Ao = 3.9, e = 0.1. 

In Fig. |l| we show the space-time picture of the active 
sites for four different sets of parameters. In all cases 
we take g{X) uniform within the range (Aq, Ao + A A) and 
zero outside, with A A = 0.1. This four cases can be 
considered paradigmatic of the behaviours observed in 
a wide region of parameter space. Fig. 0(a) and (b) 



correspond to a rather high value of the coupling, e = 0.5. 
Fig. |l|(c) and (d) have, instead, e = 0.1. The value 
of Ao is 3.9 in (b) and (d), 3.7 in (a), and 3.8 in (c). 
Depending on the value of the parameters, the behaviour 
of the successive minima can be classified in one of the 
following four categories. 

(a) discontinuous lines at i = constant: The activity 
is concentrated in a few sites of the system, and in 
first approximation the same sites remain active as 
time goes by. 

(b) uniform: The activity is spread all over the system 
without any apparent order. 

(c) "worms" : The activity is mostly concentrated in a 
few sites of the system. The activity, however, wan- 
ders as time goes by following a random pattern. 

(d) clustcrized: The activity is spread all over the sys- 
tem, but in clusters. At variance with (b), at any 
given time, one sees regions in space where there is 
no activity. 

Let us discuss in more detail Fig. gc). The worms are 
created during the transient. Once created, each worm 
wanders in space till it encounters a second worm. At 
that moment they annihilate each other. We have ob- 
served that the inactive regions correspond to a periodic 
pattern in space psf . The worms, in turn, correspond to 
defects in the periodicity. Eventually, the last two surviv- 
ing worms merge, and from there on the activity remains 
concentrated on a single "fat" worm. 

For each one of these cases, we have computed the first 
return time distribution. This corresponds to the distri- 
bution of times between consecutive activity in the same 
site. We observe that when the successive minima are 
distributed uniformly in space, the first return time dis- 
tribution is close to an exponential (see Fig. ^(b)). As 
soon as the activity is not uniformly distributed, be it 
lines or irregular clusters, the first return time distribu- 
tion exhibits a power law decay (asymptotically in case 
(d)). 

To make this picture more quantitative, we present 
in Fig. ^ a phase diagram, in A-e space (AA is held 
constant). Phase I is characterised by the presence of 
power law behaviour in the first return time distribu- 
tion, P{t) ^ T~". The value of the exponent of this 
distribution is not constant within the region. Phase II, 
on the other hand, is characterised by a non power law 
behaviour (close to exponential) in the first return time 
distribution (see Fig. |(b)). This is tantamount to saying 
that in phase I there is clusterization of the activity, that 
is absent in phase II. In Fig. |^ we indicate also a third 
phase (III), which corresponds to an asymptotic power 
law (see Fig. ^(d)), with an exponent that is different 
from that in region I. The border of this region (that 
we show only schematically) is not as well defined as the 
border between regions I and II. 
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FIG. 2. Distribution of first return time. Tlie parame- 
ters correspond to tiiose in Fig. |f|: TV = 500, AA = O.I, 
(a)Ao = 3.7, e — 0.5. The dashed line corresponds to a power 
law with exponent a = 1.80(3). (b)Ao = 3.9, e = 0.5. The 
dashed line, shown as a reference, corresponds to an exponen- 
tial decay, (c) Ao = 3.8, e — O.I. The dashed line corresponds 
to a power law with exponent a — 1.57(2). (d)Ao = 3.9, 
e = O.I. The dashed line corresponds to a power law with 
exponent a — 1.53(1), fitted to the asymptotic region. 

This picture holds, qualitatively, even in the absence of 
the evolutionary dynamics. Indeed, if in every time step 
we track the position of the minimum, but do not change 
the value of its A (this corresponds, effectively, to setting 
A A = 0), we still obtain the four abovementioned regimes 
and the corresponding first return time distributions. As 
AA changes, the value of the exponent in the first return 
time distribution (in region I) changes as well. Moreover, 
the area of region II increases with increasing AA. 

The presence of evolutionary dynamics (i. e. AA > 0), 
has yet another effect. The distribution of A evolves 
in time, from an originally uniform to a stationary 
non-uniform one. Indeed, this stationary distribution 
is peaked close to Aq and monotonically decresases for 
larger values of A: The system has "self-organized". Since 
the extremal dynamics favors higher values of x, and the 
higher the value of A the more likely small values of x 
are Q , then large values of A are more likely to be up- 
dated. The exact shape of this distribution depends on 
the values of the parameters . 

Summarizing, the dynamics of extremal coupled maps 
exhibits a behavior usually associated with criticality: 

• First return time distribution is a power law (even 
in the absence of extremal dynamics) 



• there is self-organization in the A space 

• the active regions are localized in space, very much 
like the avalanches in SOC models (BS). 
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FIG. 3. Phase diagram in A-e space. iV = 200, AA = O.I. 
Averaged over 10 realizations. 

It is worth emphasizing that this behavior (except the 
self organization in A space) remains even in the absence 
of extremal dynamics. In particular, the first return time 
distribution exhibits a power law behavior. This implies 
that power law distributions are compatible with non- 
extremal systems, as has been previously observed by 
Newman an coworkers , in the context of noise-driven 
systems. 

Both SOC and CML models have been used sepa- 
rately in the past to describe evolving ecosystems. To our 
knowledge, the model presented here is the first synthe- 
sis of both ideas that includes the level structure (geno- 
type/phenotype) of real organisms. 



* E-mail: abramson@cab. cnea.gov. ar. Present address: 
Centro Atomico Bariloche, 8400 Bariloche, Argentina. 

^ E-mail: jose@mpipks-dresden.mpg.de. 
[I] J.M. Carlson and J.S. Langer, Phys. Rev. Lett. 62, 2632 
(1989). 

[2] L. Knopoffand D. Sornette, J. Physique I 5, 1682 (1995). 

[3] S.I. Zaitsev, Physica A 189, 411 (1992). 

[4] A. Petri, G. Paparo, A. Vespignani, A. Alippi and M. 

Costantini, Phys. Rev. Lett. 73, 3423 (1994). 
[5] P. Diodati, F. Marchesoni and S. Piazza, Phys. Rev. Lett. 

67, 2239 (I99I). 
[6] G. CaldareUi, F. di ToUa and A. Petri, Phys. Rev. Lett. 

77, 2503 (1996). 
[7] D. Wilkinson and J.F. Willemsen, J. Phys. A 16, 3365 

(1983). 

[8] M. Cieplak and M.O. Robbins, Phys. Rev. Lett. 60, 2042 
(1988). 

[9] K. Sneppen, Phys. Rev. Lett. 69, 3539 (1992). 

[10] K. Sneppen, Phys. Rev. Lett. 71, lOI (1993). 

[II] A. Rinaldo, I. Rodriguez-Iturbe, R. Rigon, E. Ijjasz- 
Vasquez and R.L. Bras, Phys. Rev. Lett. 70, 822 (1993). 

[12] A. Maritan, F. Colaiori, A. Flammini, M. Cieplak and 
J.R. Banavar, Science 272, 984 (1996). 

[13] G. CaldareUi, A. Giacometti, A. Maritan, I. Rodriguez- 
Iturbe and A. Rinaldo, Phys. Rev. E 55, R4865 (1997). 



3 



[14] M.D. Raup, Science 251, 1530 (1986). 

[15] N. Eldredge and S.J. Gould, Punctuated equilibna: an 
alternative to phyletic gradualism in Models in Paleobiol- 
ogy, Freeman, Cooper and Co., San Francisco (1972). 

[16] S.J. Gould and N. Eldredge, Paleobiology 3, 115 (1977). 

[17] N. Eldredge and S.J. Gould, Nature 332, 211 (1988). 

[18] P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. Lett. 59, 
381 (1987). 

[19] P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. A 38, 
364 (1988). 

[20] M. Paczuski, S. Maslov and P. Bak, Phys. Rev. E 53, 
414 (1996). 

[21] P. Bak and K. Sneppen, Phys. Rev. Lett. 71, 4083 (1993). 
[22] P. De Los Rios, A. Valleriani and J. L. Vega, Phys. Rev. 

E 56, 4876 (1997). 
[23] P. De Los Rios, A. Valleriani and J. L. Vega, Phys. Rev. 

E 57, 6451 (1998). 
[24] S. Maslov, M. Paczuski and P. Bak, Phys. Rev. Lett. 73, 

2162 (1994). 

[25] P. Grassberger, Phys. Lett. A 200, 277 (1995). 

[26] The logistic map is a paradigmatic model of the evolution 
of Xi in the context of population dynamics. As such, it 
has been carefully studied and many its properties are 
well known. 

[27] Coupled Map Lattices have been extensively studied by 
Kaneko and coworkers , as well as many other authors 
(see and references therein) . They exhibit frozen pat- 
terns as well as travelling waves depending on the values 
of e and g{X). 

[28] K. Kaneko, Theory and applications of coupled map lat- 
tices (John Wiley & Sons, Chichester, 1993). 

[29] H. Chate and P. Manneville, Prog. Theor. Phys. 87, 1 
(1992). 

[30] H. G. Schuster, Determistic Chaos VCH Verlag, Wein- 

heim, 1989). 
[31] This topic is currently under investigation. 
[32] B. W. Roberts and M. E. J. Newman, J. Theor. Biol. 

180, 39 (1996). 



4 



